Ham Radio 2000 #1/Ham Radio 2000.iso/ham2000/misc/worldmap/mapvu20/mapproj.pas

Pascal/Delphi Source File | 1989-01-11 | 32.3 KB | 829 lines

Unit mapproj;  { projection stuff for MapView }

{ Copyright A.J. van den Bogert and Gisbert W.Selke Dec 1988                 }

{$UNDEF  DEBUG }                       { DEFINE while debugging }

{$IFDEF DEBUG }
{$A+,R+,S+,I+,D+,F-,V-,B-,L+ }      { turn checking on while debugging }
{$ELSE }
{$A+,R-,S-,I+,D-,F-,V-,B-,L+ }
{$ENDIF }
{$M 65500,65500,560000}

{$IFDEF CPU87 }
{$N+ }
{$ELSE }
{$N- }
{$ENDIF }

Interface

Uses  Crt, Graph, mapgraph;

Const maxbuf        = 100;      { max. size of curve buffer }
      trigtablesize = 9000;     { size of table for precomputed trig values }
      epsilon       = 1e-5;     { nearly 0.0 - for comparison of reals }

Type buffer     = Array [1..maxbuf] Of real;
     trigtable  = Array [0..trigtablesize] Of real;
     contbuffer = Array [1..maxbuf] Of boolean;
     ptypes = (none, mercator, ortho, lambert, azinorth, azisouth);

Var projtype : ptypes;              { projection type }
    fasttrig : boolean;             { do we have trig tables? }
    midlon, midlat : real;          { median point for ortho/lambert }
    lonmin, lonmax, latmin, latmax : real; { window borders }


Function  rmin(x, y : real) : real;          { minimum of two reals }
Function  rmax(x, y : real) : real;          { maximum of two relas }
Function  isignum(i : integer) : integer;    { sign of input (integer) }
Function  rsignum(x : real) : integer;       { sign of input (real) }
Function  scalex(x : real) : integer;        { scale in x direction }
Function  scaley(y : real) : integer;        { scale in y direction }
Function  mercproj(lat:real) : real;         { Mercator projection }
Function  invmercproj(y:real) : real;        { inverse of same }
Procedure orthoproj(lon,lat:real; Var xg,yg,zg:real);{ orthogr. projection }
Procedure orthorot(xg,yg,zg:real; Var x,y,z:real);   { globe rotation }
Procedure lambproj(lon,merclat:real; Var x,y:real);  { Lambert projection }
Procedure aziproj(lon,lat:real; Var x,y:real);       { azimuthal projection }
Procedure project(lon,lat:real; Var x,y:real);       { general projection }
Procedure invproject(x,y:real; Var lon,lat:real);    { inverse of them all }
Procedure checkwindow;                       { correct window, if necessary }
Procedure adaptscale;                        { expand limits to fill screen }
Procedure getpoint(Var x,y:real;             { get coord. via crosshairs }
                   Var finish,firstin:boolean; endgetpoint:boolean);
Procedure drawgrid(gridlon,gridlat:real; grcolour:integer;
                   smile:boolean);           { draw grid on globe }
Procedure mapborder(brcolour:integer; smile:boolean);{ determine&draw borders }
Procedure setprojtype(pt:ptypes);            { set projection parameters }
Procedure readtrigs(trigname : string);      { try to read trig tables }
Procedure dispotrigs;                        { dispose of trig tables }


Implementation

Var
 xlo, xhi, xppu, ylo, yhi, yppu : real;     { scaling parameters }
 zproj : real;                              { z-coordinate of transformation }
 raddeg, pihalf, pifourth, pidegree : real; { Pi-related constants }
 p1,p2,p3,p4,p5,p6,p7,p8,p9 : real;         { projection matrix }
 q1,q2,q3,q4,q5,q6,q7,q8,q9 : real;         { inverse proj. matrix }
 mercptr, sinptr : ^trigtable;            { precomputed Mercator / sine values }

Function rmin(x, y : real) : real;
{ minimum of arguments                                                       }
Begin                                                                 { rmin }
  If x < y Then rmin := x Else rmin := y;
End;                                                                  { rmin }

Function rmax(x, y : real) : real;
{ maximum of arguments                                                       }
Begin                                                                 { rmax }
  If x > y Then rmax := x Else rmax := y;
End;                                                                  { rmax }

Procedure hilo(Var x, y : buffer; Var c : contbuffer; n : integer);
{ update x and y limits of curve buffer                                      }
  Var i : integer;
Begin                                                                 { hilo }
  For i := 1 to n Do
  Begin
    If c[i] Then
    Begin
      If x[i] < xlo Then xlo := x[i];
      If x[i] > xhi Then xhi := x[i];
      If y[i] < ylo Then ylo := y[i];
      If y[i] > yhi Then yhi := y[i];
    End;
  End;
End;                                                                  { hilo }

Function isignum(i : integer) : integer;
{ -1, if negative; +1, if positive; 0, if 0                                  }
Begin                                                              { isignum }
  If i = 0 Then isignum := 0
           Else If i > 0 Then isignum := 1 Else isignum := -1;
End;                                                               { isignum }

Function rsignum(x : real) : integer;
{ -1, if negative; +1, if positive; 0, if 0                                  }
Begin                                                              { rsignum }
  If Abs(x) < epsilon Then rsignum := 0
                      Else If x > 0.0 Then rsignum := 1 Else rsignum := -1;
End;                                                               { rsignum }

Function scalex(x : real) : integer;
{ scale x value to screen dimension                                          }
Begin                                                               { scalex }
  scalex := Trunc(xppu * (x - xlo));
End;                                                                { scalex }

Function scaley(y : real) : integer;
{ scale y value to screen dimension                                          }
Begin                                                               { scaley }
  scaley := Trunc(yppu * (yhi - y));
End;                                                                { scaley }

Function tan(x : real) : real;
{ Niklaus forgot it, tsk, tsk                                                }
Begin                                                                  { tan }
  tan := Sin(x) / Cos(x);
End;                                                                   { tan }

Function  arcsin(x : real) : real;
Begin                                                               { arcsin }
  arcsin := ArcTan(x / Sqrt(1.0 - Sqr(x)));
End;                                                                { arcsin }

Function fastlntan(x : real) : real;
{ get Ln(Tan(..)) out of precomputed table                                   }
Begin                                                            { fastlntan }
  If x >= 0.0 Then fastlntan :=  mercptr^[Round(x*100.0)]
              Else fastlntan := -mercptr^[Round(-x*100.0)];
End;                                                             { fastlntan }

Function fastsin(x : real) : real;
{ get Sin(x) out of precomputed table                                        }
  Var ix : integer;
Begin                                                              { fastsin }
  While x >=  180.0 Do x := x - 360.0;
  While x <= -180.0 Do x := x + 360.0;
  ix := Round(x*100.0);
  If ix >= 0 Then If ix <=  9000 Then fastsin :=  sinptr^[      ix]
                                 Else fastsin :=  sinptr^[18000-ix]
             Else If ix >= -9000 Then fastsin := -sinptr^[     -ix]
                                 Else fastsin := -sinptr^[18000+ix];
End;                                                               { fastsin }

Function fastcos(x : real) : real;
{ get Cos(x) out of precomputed table                                        }
  Var ix : integer;
Begin                                                              { fastcos }
  While x >=  180.0 Do x := x - 360.0;
  While x <= -180.0 Do x := x + 360.0;
  ix := Round(x*100.0);
  If ix >= 0 Then If ix <=  9000 Then fastcos :=  sinptr^[ 9000-ix]
                                 Else fastcos := -sinptr^[-9000+ix]
             Else If ix >= -9000 Then fastcos :=  sinptr^[ 9000+ix]
                                 Else fastcos := -sinptr^[-9000-ix];
End;                                                               { fastcos }

Function mercproj(lat : real) : real;
{ compute Mercator projection of latitude                                    }
Begin                                                             { mercproj }
  If fasttrig Then mercproj := fastlntan(lat)
  Else
  Begin
    { leave it alone if too close to pole }
    If (Abs(lat) > 85.0) Then mercproj := lat
                         Else mercproj := Ln(Tan(pidegree*lat + pifourth));
  End;
End;                                                              { mercproj }

Function invmercproj(y : real) : real;
{ compute inverse of Mercator projection of latitude                         }
Begin                                                          { invmercproj }
  If (Abs(y) > 85.0) Then invmercproj := y
                     Else invmercproj := (ArcTan(Exp(y)) - pifourth)/pidegree;
End;                                                           { invmercproj }

Procedure orthoproj(lon, lat : real; Var xg, yg, zg : real);
{ calculate R3 coordinates of given point                                    }
  Var coslat : real;
Begin                                                            { orthoproj }
  If fasttrig Then
  Begin
    coslat := fastcos(lat);
    xg := fastcos(lon)*coslat;
    yg := fastsin(lon)*coslat;
    zg := fastsin(lat);
  End
  Else
  Begin
    lon := lon*raddeg;
    lat := lat*raddeg;
    coslat := Cos(lat);
    xg := Cos(lon)*coslat;
    yg := Sin(lon)*coslat;
    zg := Sin(lat);
  End;
End;                                                             { orthoproj }

Procedure invorthoproj(xg, yg, zg : real; Var lon, lat : real);
{ calculate longitude and latitude                                           }
  Var r : real;
Begin                                                         { invorthoproj }
  r := Sqr(xg) + Sqr(yg);
  r := SqRt(rmax(Sqr(xg)+Sqr(yg),0.0));
  If r > epsilon Then lat := ArcTan(zg/r)/raddeg
                 Else lat := rsignum(zg) * 90.0;
  If Abs(xg) > epsilon Then
  Begin
    lon :=  ArcTan(yg/xg)/raddeg;
    If xg < 0.0 Then If yg >= 0.0 Then lon := lon + 180.0
                                  Else lon := lon - 180.0;
  End Else lon := rsignum(yg) * 90.0;
End;                                                          { invorthoproj }

Procedure orthorot(xg, yg, zg : real; Var x, y, z : real);
{ rotate globe                                                               }
Begin                                                             { orthorot }
  x := p4*xg + p5*yg;
  y := p7*xg + p8*yg + p9*zg;
  z := p1*xg + p2*yg + p3*zg;
End;                                                              { orthorot }

Procedure lambproj(lon, merclat : real; Var x, y : real);
{ calculate conformal Lambert mapping of point                               }
  Var alpha, dist : real;
Begin                                                             { lambproj }
  dist := p7 * Exp(-p6 * merclat);
  If fasttrig Then
  Begin
    alpha :=    p6 * (lon-midlon);
    x     :=  dist * fastsin(alpha);
    y     := -dist * fastcos(alpha);
  End
  Else
  Begin
    alpha :=    p4 * (lon-midlon);
    x     :=  dist * Sin(alpha);
    y     := -dist * Cos(alpha);
  End;
End;                                                              { lambproj }

Procedure invlambproj(x, y : real; Var lon, lat : real);
{ calculate inverse of conformal Lambert mapping                             }
  Var alpha, dist : real;
Begin                                                          { invlambproj }
  If (Abs(x) <= epsilon) And (Abs(y) <= epsilon) Then
  Begin
    lat := 90.0;
    lon :=  0.0;
  End Else
  Begin
    If Abs(y) >= epsilon Then
    Begin
      alpha := ArcTan(-x/y);
      dist  := -y / Cos(alpha);
    End Else
    Begin
      dist := Abs(x);
      alpha := rsignum(x) * 90.0;
    End;
    lon := alpha/p4 + midlon;
    lat := (pifourth - ArcTan(Exp(Ln(dist/p7)/p6)))/pidegree;
  End;
End;                                                           { invlambproj }

Procedure aziproj(lon, lat : real; Var x, y : real);
{ calculate coordinates under area-prserving azimuthal projection            }
  Var dist : real;
Begin                                                              { aziproj }
  If fasttrig Then
  Begin
    dist := 2.0 * fastsin(45.0 - 0.5*lat);
    x    := dist * fastsin(lon);
    y    := -dist * fastcos(lon);
  End Else
  Begin
    dist := 2.0 * Sin(pifourth - lat*pidegree);
    lon  := lon * raddeg;
    x    := dist * Sin(lon);
    y    := -dist * Cos(lon);
  End;
End;                                                               { aziproj }

Procedure invaziproj(x, y : real; Var lon, lat : real);
{ inverse azimuthal projection                                               }
  Var dist : real;
Begin                                                           { invaziproj }
  If Abs(y) > epsilon Then
  Begin
    lon  := ArcTan(-x/y) / raddeg;
    If y > 0.0 Then If x < 0.0 Then lon := lon - 180.0
                               Else lon := lon + 180.0;
    dist := SqRt(Sqr(x) + Sqr(y));
  End Else
  Begin
    lon  := rsignum(x) * 90.0;
    dist := Abs(x);
  End;
  lat := 90.0 - ArcSin(0.5*dist) / pidegree;
End;                                                            { invaziproj }

Procedure project(lon, lat : real; Var x, y : real);
{ convert longitude and latitude to rectangular coordinates                  }
  Var xg, yg, zg: real;   { coordinates of point on globe }
Begin                                                              { project }
  Case projtype Of
    none:      Begin
                 x := lon;
                 y := lat;
               End;
    mercator:  Begin
                 x := lon;
                 y := mercproj(lat);
               End;
    ortho:     Begin
                 orthoproj(lon,lat,xg,yg,zg);
                 { rotate the globe }
                 orthorot(xg,yg,zg,x,y,zproj);
               End;
    lambert:   lambproj(lon,mercproj(lat),x,y);
    azinorth:  aziproj(lon,lat,x,y);
    azisouth:  aziproj(180.0-lon,-lat,x,y);
  End; { Case }
End;                                                               { project }

Procedure invproject(x, y : real; Var lon, lat : real);
{ inverse projection                                                         }
  Var xg, yg, zg : real;            { coordinates on globe }
Begin                                                           { invproject }
  Case projtype Of
    none:      Begin
                 lon := x;
                 lat := y;
               End;
    mercator:  Begin
                 lon := x;
                 lat := invmercproj(y);
               End;
    ortho:     Begin
                 zproj := SqRt(rmax(1.0-Sqr(x)-Sqr(y),0.0));
                 { rotate the globe }
                 xg := q1*zproj + q2*x + q3*y;
                 yg := q4*zproj + q5*x + q6*y;
                 zg := q7*zproj +        q9*y;
                 invorthoproj(xg,yg,zg,lon,lat);
               End;
    lambert:   invlambproj(x,y,lon,lat);
    azinorth:  invaziproj(x,y,lon,lat);
    azisouth:  Begin
                 invaziproj(x,y,lon,lat);
                 If lon >= 0.0 Then lon :=  180.0 - lon
                               Else lon := -180.0 - lon;
                 lat := -lat;
               End;
  End; {case}
End;                                                            { invproject }

Procedure checkwindow;
{ adjust window to avoid log(0) and other nasty things                       }

  Procedure rswap(Var a, b : real);
  { swap two reals                                                           }
    Var temp : real;
  Begin                                                              { rswap }
    temp := a;
    a := b;
    b := temp;
  End;                                                               { rswap }

Begin                                                          { checkwindow }
  If latmin > latmax Then rswap(latmin,latmax);
  If lonmin > lonmax Then rswap(lonmin,lonmax);
  If Abs(latmax-latmin) < epsilon Then
  Begin
    latmax := latmax + 1.0;
    latmin := latmin - 1.0;
  End;
  If Abs(lonmax-lonmin) < epsilon Then
  Begin
    lonmax := lonmax + 1.0;
    lonmin := lonmin - 1.0;
  End;
End;                                                           { checkwindow }

Procedure adaptscale;
{ expand one of the directions to fill entire screen                         }
  Var mean, scalfac, mlatmin, mlatmax : real;
Begin                                                           { adaptscale }
  Case projtype Of
    none : Begin
             scalfac := xmaxpix / (ymaxpix*aspect);
             If (lonmax-lonmin)/(latmax-latmin) > scalfac Then
             Begin
               { X scale larger than Y: increase Y }
               mean := (latmin + latmax) * 0.5;
               latmin := mean - (lonmax - lonmin) * 0.5 / scalfac;
               latmax := mean + (lonmax - lonmin) * 0.5 / scalfac;
             End Else
             Begin
               { Y scale larger than X: increase X }
               mean := (lonmin + lonmax) * 0.5;
               lonmin := mean - (latmax - latmin) * 0.5 * scalfac;
               lonmax := mean + (latmax - latmin) * 0.5 * scalfac;
             End;
           End;
    mercator : Begin
             checkwindow;
             scalfac := xmaxpix / (ymaxpix*aspect*raddeg);
             mlatmax := mercproj(latmax);
             mlatmin := mercproj(latmin);
             If (lonmax-lonmin)/(mlatmax-mlatmin) > scalfac Then
             Begin
               { X scale larger than Y: increase Y }
               mean := (mlatmin + mlatmax) * 0.5;
               latmin := invmercproj(mean - (lonmax-lonmin)*0.5/scalfac);
               latmax := invmercproj(mean + (lonmax-lonmin)*0.5/scalfac);
               checkwindow;
             End Else
             Begin
               { Y scale larger than X: increase X }
               mean := (lonmin + lonmax) * 0.5;
               lonmin := mean - (mlatmax-mlatmin)*0.5*scalfac;
               lonmax := mean + (mlatmax-mlatmin)*0.5*scalfac;
             End;
           End;
    ortho    : ;  { not applicable }
    lambert  : ;  { not applicable }
    azinorth : ;  { not applicable }
    azisouth : ;  { not applicable }
  End;
End;                                                            { adaptscale }

Procedure getlonline(lon : real; Var x, y : buffer;
                                 Var c : contbuffer; Var n : integer);
{ create rect. coordinates of line of constant longitude                     }
  Var i : integer;
      lat, step, xx, yy : real;
Begin                                                           { getlonline }
  Case projtype Of
    none, mercator, lambert,
    azinorth, azisouth : Begin
                      project(lon,latmin,x[1],y[1]);
                      project(lon,latmax,x[2],y[2]);
                      c[1] := True; c[2] := True;
                      n := 2;
                    End;
    ortho         : Begin
                      step := (latmax-latmin) / Pred(maxbuf);
                      lat := latmin;
                      For i := 1 To maxbuf Do
                      Begin
                        project(lon,lat,xx,yy);
                        x[i] := xx;
                        y[i] := yy;
                        c[i] := zproj > -epsilon;
                        lat := lat + step;
                      End;
                      n := maxbuf;
                    End;
  End; {case}
End;                                                            { getlonline }

Procedure getlatline(lat : real; Var x, y : buffer;
                                 Var c : contbuffer; Var n : integer);
{ create rect. coordinates of line of constant latitude                      }
  Var i: integer;
      lon, step, xx, yy : real;
Begin                                                           { getlatline }
  Case projtype Of
    none, mercator: Begin
                      n := 2;
                      project(lonmin,lat,x[1],y[1]);
                      project(lonmax,lat,x[2],y[2]);
                      c[1] := True; c[2] := True;
                    End;
    ortho, lambert,
    azinorth, azisouth : Begin
                      zproj := 1.0; { dummy for all except ortho }
                      If Abs(lat) >= 89.0 Then
                      Begin
                        step := 360.0  / Pred(maxbuf);
                        lon  := midlon - 180.0;
                      End
                      Else Begin
                        step := (lonmax-lonmin) / Pred(maxbuf);
                        lon  := lonmin;
                      End;
                      For i := 1 to maxbuf Do
                      Begin
                        project(lon,lat,xx,yy);
                        x[i] := xx;
                        y[i] := yy;
                        c[i] := zproj > -epsilon;
                        lon := lon + step;
                      End;
                      n := maxbuf;
                    End;
  End; {case}
End;                                                            { getlatline }

Procedure getpoint(Var x, y : real;
                   Var finish, firstin : boolean; endgetpoint : boolean);
{ get a coordinate pair given by crosshairs                                  }
  Var i, j : integer;
      c : char;
Begin                                                             { getpoint }
  i := Round(xppu*(x - xlo));
  j := Round(yppu*(yhi - y));
  If firstin Then
  Begin
    vline(i);
    hline(j);
    firstin := False;
  End;
  Repeat
    c := ReadKey;
    If (c = #0) Then c := ReadKey;
    Case c Of
     uparr  : If j > 0 Then Begin hline(j); Dec(j); hline(j); End;
     dnarr  : If j < ymaxpix Then Begin hline(j); Inc(j); hline(j); End;
     lfarr  : If i > 0 Then Begin vline(i); Dec(i); vline(i); End;
     rtarr  : If i < xmaxpix Then Begin vline(i); Inc(i); vline(i); End;
     '8', cuparr : If j > 10 Then
                               Begin hline(j); j := j-10; hline(j); End;
     '2', cdnarr : If j < ymaxpix-10 Then
                               Begin hline(j); j := j+10; hline(j); End;
     '4', clfarr : If i > 10 Then
                               Begin vline(i); i := i-10; vline(i); End;
     '6' ,crtarr : If i < xmaxpix-10 Then
                               Begin vline(i); i := i+10; vline(i); End;
    End;
    finish := c In [' ','q','Q',ctrlc,cr,esc];
  Until finish Or Not Endgetpoint;
  If finish Then Begin
    vline(i); hline(j);
  End;
  x := xlo + i/xppu;
  y := yhi - j/yppu
End;                                                              { getpoint }

Procedure curve(Var x, y: buffer; Var c : contbuffer; n: integer);
{ draw curve of n points                                                     }
  Var xo, yo, xn, yn, i : integer;
Begin                                                                { curve }
  xo := scalex(x[1]); yo := scaley(y[1]);
  For i := 2 to n Do
  Begin
    xn := scalex(x[i]); yn := scaley(y[i]);
    If c[Pred(i)] And c[i] Then Line(xo,yo,xn,yn);
    xo := xn; yo := yn;
  End;
End;                                                                 { curve }

Procedure dotcurve(Var x, y: buffer; Var c : contbuffer; n:integer);
{ draw dotted curve of n points                                              }
  Var xo, yo, xn, yn, i : integer;
      dotflag : boolean;
Begin                                                             { dotcurve }
  dotflag := True;
  xo := scalex(x[1]); yo := scaley(y[1]);
  For i := 2 to n Do
  Begin
    xn := scalex(x[i]); yn := scaley(y[i]);
    If c[Pred(i)] And c[i] Then dotline(xo,yo,xn,yn,dotflag);
    xo := xn; yo := yn;
  End;
End;                                                              { dotcurve }

Procedure drawgrid(gridlon, gridlat : real; grcolour : integer;
                   smile : boolean);
{ display grid according to set intervals                                    }
  Var lon, lat : real;
      x, y : buffer;
      c : contbuffer;
      n : integer;
      dotfl : boolean;
Begin                                                             { drawgrid }
  SetColor(word(Abs(grcolour)));
  dotfl := grcolour < 0;
  lon := gridlon * Trunc(lonmin/gridlon);
  If lon < lonmin Then lon := lon + gridlon;
  While lon < lonmax Do
  Begin
    If smile Then showprogress(1);
    getlonline(lon,x,y,c,n);
    If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
    lon := lon + gridlon;
  End;
  lat := gridlat * Trunc(latmin/gridlat);
  If lat < latmin Then lat := lat + gridlat;
  While lat < latmax Do
  Begin
    If smile Then showprogress(1);
    getlatline(lat,x,y,c,n);
    If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
    lat := lat + gridlat;
  End;
  SetColor(colourglb);
End;                                                              { drawgrid }

Procedure mapborder(brcolour : integer; smile : boolean);
{ determine scale and draw border of map area                                }

  Var x, y, xa, ya  : buffer;
      c, ca : contbuffer;
      n, na : integer;

  Procedure mapextremes;
  { determine x and y extremes                                               }
    Var lon, lat, x1, y1, rfact : real;
        i : integer;
  Begin                                                        { mapextremes }
    xlo := 1e30;
    xhi := -1e30;
    ylo := 1e30;
    yhi := -1e30;
    Case projtype Of
      none, mercator, lambert,
      azinorth, azisouth : Begin
          getlatline(latmin,x,y,c,n);
          getlatline(latmax,xa,ya,ca,na);
        End;
      ortho : Begin
          getlonline(lonmin,x,y,c,n);
          getlonline(lonmax,xa,ya,ca,na);
          hilo(x,y,c,n);
          hilo(xa,ya,ca,na);
          If smile Then showprogress(1);
          getlatline(latmin,x,y,c,n);
          getlatline(latmax,xa,ya,ca,na);
          If smile Then showprogress(1);
          hilo(x,y,c,n);
          hilo(xa,ya,ca,na);
          { calculate horizon }
          rfact := 2.0 / Pred(maxbuf);
          x1 := -1.0;
          For i := 1 To maxbuf Do
          Begin
            y1 := 1.0 - Sqr(x1);
            If y1 > 0.0 Then y1 := SqRt(y1) Else y1 := 0.0;
            invproject(x1,y1,lon,lat);
            x[i] := x1; y[i] := y1;
            c[i] :=  (lon >= lonmin-epsilon) And (lon <= lonmax+epsilon) And
                     (lat >= latmin-epsilon) And (lat <= latmax+epsilon);
            invproject(x1,-y1,lon,lat);
            xa[i] := x1; ya[i] := -y1;
            ca[i] := (lon >= lonmin-epsilon) And (lon <= lonmax+epsilon) And
                     (lat >= latmin-epsilon) And (lat <= latmax+epsilon);
            x1 := x1 + rfact;
          End;
          n := maxbuf; na := maxbuf;
          If smile Then showprogress(1);
        End;
    End;
    hilo(x,y,c,n);
    hilo(xa,ya,ca,na);
    If xhi - xlo < epsilon Then xhi := xhi + 1.0;
    If yhi - ylo < epsilon Then yhi := yhi + 1.0;
    If smile Then showprogress(1);
  End;                                                         { mapextremes }

  Procedure drawborder;
  { draw border                                                              }
    Var dotfl : boolean;
  Begin                                                         { drawborder }
    SetColor(word(Abs(brcolour)));
    dotfl := brcolour < 0;
    If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
    If dotfl Then dotcurve(xa,ya,ca,na) Else curve(xa,ya,ca,na);
    If smile Then showprogress(1);
    If projtype In [ortho,azinorth,azisouth] Then
    Begin
      If lonmax - lonmin < 360.0 - epsilon Then
      Begin
        getlonline(lonmin,x,y,c,n);
        If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
        getlonline(lonmax,x,y,c,n);
        If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
      End;
      If smile Then showprogress(1);
      getlatline(latmin,x,y,c,n);
      getlatline(latmax,xa,ya,ca,na);
      If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
      If dotfl Then dotcurve(xa,ya,ca,na) Else curve(xa,ya,ca,na);
    End
    Else Begin
      getlonline(lonmin,x,y,c,n);
      If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
      getlonline(lonmax,x,y,c,n);
      If dotfl Then dotcurve(x,y,c,n) Else curve(x,y,c,n);
    End;
    SetColor(colourglb);
  End;                                                          { drawborder }

Begin                                                            { mapborder }
  mapextremes;
  { set scale parameters (pixels per unit), ensure reasonable aspect ratio }
  xppu := xmaxpix / (xhi-xlo);
  yppu := ymaxpix / (yhi-ylo);
  Case projtype Of
    none, ortho, lambert,
    azinorth, azisouth :  { equal scale in x and y }
             Begin
               xppu := rmin(xppu,yppu*aspect);
               yppu := rmin(yppu,xppu/aspect);
             End;
    mercator : { ppu ratio y/x must be pi/180 }
             Begin
               xppu := rmin(xppu,yppu*aspect*raddeg);
               yppu := rmin(yppu,xppu/(aspect*raddeg));
             End;
  End; {case}
  drawborder;
End;                                                             { mapborder }

Procedure setprojtype(pt : ptypes);
{ set projection type; calculate parameters                                  }

  Procedure lambmatrix;
  { calculate matrix elements for Lambert projection and its inverse         }
  Begin                                                         { lambmatrix }
    p7 := pihalf - midlat*raddeg;
    p6 := Cos(p7);
    p4 := p6 * raddeg;
    p7 := Tan(p7) / Exp(p6 * Ln(Tan(p7*0.5)));
  End;                                                          { lambmatrix }

  Procedure orthomatrix;
  { calculate matrix elements for orthogonal projection and its inverse      }
  Begin                                                        { orthomatrix }
    p1 := Cos(midlon*raddeg)*Cos(midlat*raddeg);
    p2 := Sin(midlon*raddeg)*Cos(midlat*raddeg);
    p3 := Sin(midlat*raddeg);
    p4 := -Sin(midlon*raddeg);
    p5 := Cos(midlon*raddeg);
    p6 := 0.0;
    p7 := -Cos(midlon*raddeg)*Sin(midlat*raddeg);
    p8 := -Sin(midlon*raddeg)*Sin(midlat*raddeg);
    p9 := Cos(midlat*raddeg);
    { inverse matrix is transpose}
    q1 := p1;
    q2 := p4;
    q3 := p7;
    q4 := p2;
    q5 := p5;
    q6 := p8;
    q7 := p3;
    q8 := 0.0;
    q9 := p9;
  End;                                                         { orthomatrix }

Begin                                                             { projtype }
  projtype := pt;
  Case projtype Of
    none     : ;
    mercator : ;
    ortho    : orthomatrix;
    lambert  : lambmatrix;
    azinorth, azisouth : Begin
                 If latmax >= -latmin Then projtype := azinorth
                                      Else projtype := azisouth;
                 midlon := 0.0; midlat := 0.0;
               End;
  End;
End;                                                              { projtype }

Procedure readtrigs(trigname : string);
{ read trig tables into memory, if possible                                  }
  Var trigf : File Of trigtable;
Begin                                                            { readtrigs }
  New(mercptr);
  New(sinptr);
  fasttrig := (mercptr <> Nil) And (sinptr <> Nil);
  If (mercptr <> Nil) And (sinptr = Nil) Then Dispose(mercptr);
  If fasttrig Then
  Begin
    Assign(trigf,trigname);
    {$I- } Reset(trigf);
    read(trigf,mercptr^);
    read(trigf,sinptr^);
    Close(trigf);
    {$I+ }
    fasttrig := IOResult = 0;
  End;
End;                                                             { readtrigs }

Procedure dispotrigs;
{ dispose of trig tables in memory                                           }
Begin                                                           { dispotrigs }
  If fasttrig Then
  Begin
    Dispose(mercptr);
    Dispose(sinptr);
    fasttrig := False;
  End;
End;                                                            { dispotrigs }

Begin
  fasttrig := False;
  raddeg   := Pi / 180.0;
  pihalf   := Pi / 2.0;
  pifourth := Pi / 4.0;
  pidegree := Pi / 360.0;
  projtype := mercator;
  midlon   := 0.0;
  midlat   := 0.0;
  lonmin   := -200.0;
  lonmax   :=  200.0;
  latmin   :=  -80.0;
  latmax   :=   80.0;
  xppu     := 1.0;
  yppu     := 1.0;
End.